#ifndef GUARD_Newton_Fractal
#define GUARD_Newton_Fractal

// Newton_Fractal.h
// Licensed under GPLv2
// Terms available at http://www.gnu.org/licenses/gpl-2.0.txt

#include <vector>
#include <complex>
#include "CImg.h"
#include "Complex_Polynomial.h"
#include <stdexcept>

using namespace std;
using namespace cimg_library;

complex<double> newton ( complex<double> pt, Complex_Polynomial& fn );
CImg<unsigned char> plot ( Complex_Polynomial& fn , int imgx , int imgy,
                           double xmin , double xmax, double ymin,
                           double ymax , double delta );


//First matter is to set up the actual Newton's method. This is fairly
//straight-forward.
complex<double> newton ( complex<double> pt, Complex_Polynomial& fn ) {
  //Finds the derivative, checks for division by zero, and otherwise
  //applies Newton Iteration at the current point.
  complex<double> valDeriv = fn.derivEval(pt);
  if ( abs(valDeriv) == 0 )
    throw domain_error("Divides by zero");
  return ( pt - ( fn.evalPt(pt) / valDeriv ) );
}

/*The key is to load each point by a pair of recursive for-loops and
  then iterate the newton function on it 256 times, while tracking the
  difference between the points. If the real difference or the imaginary
  difference is larger than epsilon, consider it to not converge. (By
  256 iterations it should have reached 0 anyways). Run the code in a
  try-catch block to make sure if at any point it attempts division by
  zero, the point is immediately marked black. */

CImg<unsigned char> plot ( Complex_Polynomial& fn , int imgx ,
                           int imgy, double xmin , double xmax,
                           double ymin, double ymax,
			   double delta ) {
  //Create an imgx * imgy image
  CImg<unsigned char> img (imgx , imgy , 1 , 3 );

  double xstep = ( xmax - xmin ) / imgx ;
  double ystep = ( ymax - ymin ) / imgy ;
  //  double delta = pow(10.0,-8);
  // Colours are in HSV, starting with a random number from 0 to 2pi
  // Divide circle into n even pieces, match one to each, define an
  // HSV point for each root.

  //Finding how far apart each colour is.
  int dist_around_circle = 360 / fn.roots.size() ;
  int colours [fn.roots.size()][3];
  //Randomly determined base colour.
  int base_colour = ( rand() % 359 );
  for ( vector<complex<double> >::size_type i = 0 ;
        i < fn.roots.size() ; i++ ) {
    colours[i][0] = ( base_colour + i * dist_around_circle ) % 360;
    colours[i][1] = 1;
    colours[i][2] = 1;
    // These are defined in HSV; the first is mod 360 so no problems
    // come about, while the second two are 1 so everything's max
    // saturation and value.
  }

  // This is a supplied iterator that operates over all the y values
  // in img, very handy.
  cimg_forY(img,y) {
    // Because (0,0) is in the upper-left corner, a repositioning was
    // necessary. This translates y values into y coordinates.

    double ypos = ymax - ( y * ystep );
    // An iterator on x and also some code to translate x-value to
    // x-coordinate.
    cimg_forX(img,x){

      double xpos = xmin + ( x * xstep );
      //Establishes the first point to test for iteration
      complex<double> curPt ( xpos , ypos );
      try{
        int i = 0;
        //For each point, check if the root is found
        bool rootFound = false;
        while ( i < 256 && !rootFound ) {
          for ( vector<complex<double> >::size_type j = 0 ;
                j < fn.roots.size() ; j++ ) {
            if ( abs( curPt - fn.roots[j] ) < delta ) {
              //Paint it the colour of that root

              //this uses a poorly-documented part of the code.
              img.fillC(x,y,0,colours[j][0], colours[j][1] , colours[j][2] );
              //There was a reference to that line of code somewhere
              //in the CImg docs, but I forget where.

              rootFound = true;

              //Short-circuits paired with the immediate next line
              //after the for-block
              break;
            }
          }
          //Runs the Newton iteration.
          curPt = newton( curPt, fn );

          //I used a while loop so I could check against the root being found;
          //a for-loop would have done too.
          i++;
        }
      } catch ( domain_error e ) {
        //This case applies when division by zero is attempted in the
        //derivative and shouldn't show up too often. In fact, nothing need
        //happen, because it'll just stay a black pixel. Which is my plan.
      }

    }
  }
  return img;
}

#endif
